clear all
path(path,'../TOOLS_01');

%%%%% le malha
aux=load_gmsh4('malha_32_ele.msh',[1,2]);

n_nos=aux.nbNod;
n_elementos=aux.nbTriangles;
coords=aux.POS;
topologia_aux= aux.ELE_NODES(:,1:3);
n_eletrodos=max(aux.LINES(:,3));
topologia=topologia_aux(n_eletrodos+aux.nbLines+1:end,:);

ele_tags=aux.ELE_TAGS(n_eletrodos+aux.nbLines+1:end,:);
idx_objeto=find(ele_tags(:,1)==2);
idx_malha=find(ele_tags(:,1)==1);


altura=0.05; % altura da malha
rho=2.0; %malha
rho_objeto=70.0; %objeto
rho_eletrodo=0.5;

corriente=1e-3;

rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto)=rho_objeto;
rho_malha(n_elementos+1:end)=rho_eletrodo;

topologia_eletrodo=zeros(n_eletrodos,4);
for i=1:2:n_eletrodos*2
    topologia_eletrodo((i-1)/2+1,1:3)=[aux.LINES(i,1),aux.LINES(i,2),aux.LINES(i+1,2)];
    topologia_eletrodo((i-1)/2+1,4)=topologia_aux((i-1)/2+1);
end


%%%% basal
rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto)=rho;
rho_malha(n_elementos+1:end)=rho_eletrodo;
V_basal = p_direto(corriente,n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura ,topologia_eletrodo);

%%%% objeto
rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto)=rho_objeto;
rho_malha(n_elementos+1:end)=rho_eletrodo;
V_objeto = p_direto(corriente,n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura ,topologia_eletrodo);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%  matriz de sensibilidade %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% valor inicial

rho_inicial=2.4;


rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho_inicial;
rho_malha(idx_objeto)=rho_inicial;
rho_malha(n_elementos+1:end)=rho_eletrodo;
[Y_global,Y_local] = monta_rigidez_tri(n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura);
V_basal_01 = p_direto(corriente,n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura ,topologia_eletrodo);


[ H ] = monta_jacobiano_2D(Y_global,Y_local,rho_malha,corriente,n_elementos,n_eletrodos,n_nos,coords,topologia_eletrodo);

imagem=inv(H'*H+1e-8*eye(n_elementos))*(H'*(V_objeto-V_basal_01)');


%%%%%%%%%%%%%%%%%%%% plot imagem %%%%%%%%%%%%%%%%%%%%

[ centroides ] = calc_centroides( coords, n_elementos,topologia);


max_x=max(centroides(:,1));
min_x=min(centroides(:,1));

max_y=max(centroides(:,2));
min_y=min(centroides(:,2));

res=64;
[XI,YI]=meshgrid(linspace(min_x,max_x,res),linspace(min_y,max_y,res));



F = TriScatteredInterp(centroides(:,1),centroides(:,2),imagem);
QZI=F(XI,YI);
%ZI=griddata(centroides(:,1),centroides(:,2),imagem,XI,YI);

surf(XI,YI,QZI)











